Back to Blog | LearnableLoopAI.com | Portfolio of Projects |

versioninfo() ## Julia version

Julia Version 1.10.5
Commit 6f3fdf7b362 (2024-08-27 14:19 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 12 × Intel(R) Core(TM) i7-8700B CPU @ 3.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, skylake)
Threads: 1 default, 0 interactive, 1 GC (on 12 virtual cores)

import Pkg
Pkg.add("BenchmarkTools")

using RxInfer, BenchmarkTools, Random, LinearAlgebra, Plots, LaTeXStrings, Distributions, StableRNGs

   Resolving package versions...
  No Changes to `/workspaces/Kalman filtering and smoothing/Project.toml`
  No Changes to `/workspaces/Kalman filtering and smoothing/Manifest.toml`

Pkg.status()

Status `/workspaces/Kalman filtering and smoothing/Project.toml`
  [6e4b80f9] BenchmarkTools v1.8.0
⌃ [31c24e10] Distributions v0.25.131
  [b964fa9f] LaTeXStrings v1.4.1
  [91a5bcdd] Plots v1.41.7
⌅ [86711068] RxInfer v3.10.1
⌃ [860ef19b] StableRNGs v1.0.4
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

Kalman filtering and smoothing (Part 3)

This is an analysis of the RxInfer example at https://examples.rxinfer.com/categories/basic_examples/kalman_filtering_and_smoothing/
Some symbols have been changed
Some content has been added/modified
The preference is to make the math and code names align as much as possible
Notation follows FoAI_Notation_Spec.md (v1.23)
Typeface (math)
- continuous scalars (this notebook): plain italic, e.g. $x^{\ast(t)}$, $y^{(t)}$
- continuous vectors / arrays: italic bold, e.g. $\boldsymbol{x}$, $\boldsymbol{\omega}$
- classical / engineering matrices: bold upright sans-serif, e.g. $\boldsymbol{\mathsf{A}}$, $\boldsymbol{\mathsf{\Sigma}}$
- tuples (ordered collections that may differ in shape): upright bold, e.g. $\mathbf{s}$
Code marks (align with Spec §6)
- _ prefix for a bold array/matrix when used; scalars need no _
- ː (length mark) suffix for a trajectory, e.g. xˣː, yː
- ˣ for a process/environment twin of math $ ^{} $, e.g. xˣ ↔︎ $x^{\ast}$
- ₜ for the current tick; previous tick ₜ₋₁ (Julia)
Time (math): member index in round-bracket superscripts, e.g. $x^{\ast(t)}$, $y^{(t)}$ — not a subscript
Process / environment mark: asterisk before any indices, in one superscript with time, e.g. $x^{\ast(t)}$ (KaTeX-safe; avoid the broken double-superscript form x^*_{t}^{(t)})
Named channels use subscripts; vector components use square entries $x^{(t)[n]}$ when a single array is meant

In the following set of examples the goal is to estimate hidden states of a Dynamical process where all hidden states are Gaussians.

We start our journey with a simple

1. multivariate Linear Gaussian State Space Model (LGSSM), which can be solved analytically. We then solve an
1. identification problem which does not have an analytical solution. Utimately, we show how RxInfer.jl can
1. deal with missing observations.

3 Handling Missing Data

An interesting case in filtering and smoothing problems is the processing of missing data. It can happen that sometimes your reading devices fail to acquire the data leading to missing observation.

Let us assume that the following model generates the data

\[\begin{aligned} x^{(t)} &\sim \mathcal{N}\bigl(x^{(t)}; x^{(t-1)}, w\bigr) \\ y^{(t)} &\sim \mathcal{N}\bigl(y^{(t)}; x^{(t)}, v\bigr) \end{aligned}\] where $w = 1.0$

with prior $x^{(0)} \sim \mathcal{N}(x^{(0)}; m_{x^{(0)}}, v_{x^{(0)}})$. Suppose that our measurement device fails to acquire data from time to time. In this case, instead of scalar observation $y^{(t)} \in \mathbb{R}$ we sometimes will catch missing observations.

μᵥᵥ = 0.0
σ²ᵥᵥ = 0.1

μᵥ = 0.0
σ²ᵥ = 1.0
T = 250
missing_indices = 100:125

100:125

The Generative Process

State transition function ($f_{\mathcal{E}}$)

The state transition function provides the deterministic part of the state flow. The probabilistic part is provided by the system noise:

\[\dot{x}^{\ast(t)} = f_{\mathcal{E}}(t) + w = t + w\]

where

$w \sim \mathcal{N}(w; \mu_W, \sigma^2_W)$

The *s indicate that the parameters and variables are hidden and not observed.

## state transition function
function fE(;T)
    fEː = float.(collect(1:T))
    return fEː
end
myfEː = fE(T=5)

5-element Vector{Float64}:
 1.0
 2.0
 3.0
 4.0
 5.0

Observation generation function ($g_{\mathcal{E}}$)

The observation generation function provides the deterministic part of the observation. The probabilistic part is provided by the observation noise:

\[y^{(t)} = g_{\mathcal{E}}(x^{\ast(t)}) + v = \mathrm{sin}(0.05\, x^{\ast(t)}) + v\] where

$v \sim \mathcal{N}(v; \mu_V, \sigma^2_V)$

The *s indicate that the parameters and variables are hidden and not observed.

## observation generation function
function gE(xˣː)
    gEː = sin.(0.05 * xˣː)
    return gEː
end
gE(myfEː)

5-element Vector{Float64}:
 0.04997916927067833
 0.09983341664682815
 0.14943813247359924
 0.19866933079506122
 0.24740395925452294

## Data comes from either a simulation/lab (sim|lab) OR from the field (fld)
## Data are handled either in batches (batch) OR online as individual points (point)
## Batch data accumulates either
    ## along the depth/examples dimension/axis (into the screen/page), OR
        ## typical for supervised & unsupervised learning
    ## along the time dimension/axis (down the screen page)
        ## typical for sequential decision learning (reinforcement learning & active inference)
function sim_batch_data(T; seed=123, μᵥᵥ, σ²ᵥᵥ, μᵥ, σ²ᵥ) ## simulated batch data
    rng = StableRNG(seed)
    fEː = fE(T=T)
    wː = rand(Normal(μᵥᵥ, sqrt(σ²ᵥᵥ)), T)
    xˣː = fEː + wː

    gEː = gE(xˣː)
    vː= rand(Normal(μᵥ, sqrt(σ²ᵥ)), T)
    yː = gEː + vː

    y_missingː    = similar(yː, Union{Float64, Missing}, )
    copyto!(y_missingː, yː)
    for index in missing_indices
        y_missingː[index] = missing
    end
    return fEː, xˣː, gEː, yː, y_missingː
end

sim_batch_data (generic function with 1 method)

fEː, xˣː, gEː, yː, y_missingː = sim_batch_data(
    T; 
    μᵥᵥ=μᵥᵥ, σ²ᵥᵥ=σ²ᵥᵥ, 
    μᵥ=μᵥ, σ²ᵥ=σ²ᵥ);

xˣː

250-element Vector{Float64}:
   1.481107724031414
   1.7192650416903656
   3.116297197604394
   3.5272905989350516
   4.8730989057445635
   5.9109010968117355
   6.990251532628996
   7.932995576947265
   9.44278651810198
   9.911696647642785
  11.152564962845135
  12.306470995858703
  13.066338576227425
   ⋮
 239.0145395593193
 240.2220883615564
 241.29678666667692
 242.14960946201685
 243.10148465503286
 243.83985755477784
 244.69806598692705
 245.45887089809054
 247.02989691467616
 247.94652526496418
 248.49766533263653
 250.30043535886233

yː

250-element Vector{Float64}:
 -0.5450030724167287
 -1.0033185481982012
 -0.9898949262552779
  0.3439668931501726
  0.14691926402236977
 -0.7403448214523907
  1.564502787133004
  0.7760729921271237
  0.7520463023289832
  0.24272517034321725
  1.442845080605812
 -0.40638235380508836
  1.2250965994832435
  ⋮
 -1.1727736240730384
 -0.6754135390010005
 -0.45638025817901495
  0.47088546527278163
  0.3318526979610771
 -0.4632580867204367
 -1.1946161294614976
 -1.5790623956951926
  0.24683202296990692
 -1.1848587551491903
  0.5394852136736867
 -0.1904561726450137

gEː

250-element Vector{Float64}:
  0.07398771566589066
  0.08585741769860221
  0.15518513890573996
  0.17545166393871567
  0.24125121798799845
  0.29126131654978554
  0.3424398945284677
  0.3863304046459418
  0.4547926026938593
  0.4755462042011178
  0.5291752269058843
  0.5772227589611079
  0.6078236250493826
  ⋮
 -0.5774841078568639
 -0.5271695218676266
 -0.48076880152064677
 -0.44295334204771486
 -0.39979793126555874
 -0.3656934373932446
 -0.32543082119575617
 -0.28923450776752313
 -0.21322605553431345
 -0.16824039896867374
 -0.14101575526392607
 -0.05132628410573521

tː = 1:length(yː)
p = plot(xlabel=L"t", titlefontsize=14, legendfontsize=12, guidefontsize=12,
         bottom_margin=6Plots.mm)
p = plot!(p, tː, gEː, label="True " * L"g_{\mathcal{E}}(x^{\ast(t)})", color=:blue)
p = scatter!(p, tː, yː, label=L"y^{(t)}", markersize=2, color=:orange)
plot(p)

The Generative Model

@model function smoothing(yː, x₀)
    v ~ Gamma(shape=0.001, scale=0.001)
    x_prior ~ Normal(mean=mean(x₀), var=var(x₀))
    local x
    xₜ₋₁ = x_prior
    for t in 1:length(yː)
        xː[t] ~ Normal(mean=xₜ₋₁, precision=1.0)
        yː[t] ~ Normal(mean=xː[t], precision=v)
        xₜ₋₁ = xː[t]
    end
end

constraints = @constraints begin
    q(x_prior, xː, yː, v) = q(x_prior, xː)q(v)q(yː)
end

Constraints: 
  q(x_prior, xː, yː, v) = q(x_prior, xː)q(v)q(yː)

x₀_prior = NormalMeanVariance(0.0, 1000.0)
initm = @initialization begin
    q(v) = Gamma(0.001, 0.001)
end

result = infer(
    model = smoothing(x₀=x₀_prior), 
    data  = (yː = y_missingː,), 
    constraints = constraints,
    initialization = initm, 
    returnvars = (xː = KeepLast(),),
    iterations = 20
);

plot(gEː, label="True " * L"g_{\mathcal{E}}(x^{\ast(t)})", legend=:bottomright, color=:blue,
     xlabel=L"t", titlefontsize=14, legendfontsize=12, guidefontsize=12,
     bottom_margin=6Plots.mm)
scatter!(
    missing_indices,
    gEː[missing_indices], ms=2, opacity=0.75, label="Missing region", color=:blue)
plot!(
    mean.(result.posteriors[:xː]),
    ribbon=var.(result.posteriors[:xː]), label="Estimated " * L"x^{(t)}")

Kalman filtering and smoothing (Part 3)

Kalman filtering and smoothing (Part 3)

3 Handling Missing Data

The Generative Process

State transition function (\(f_{\mathcal{E}}\))

Observation generation function (\(g_{\mathcal{E}}\))

The Generative Model